#include "stdinc.h" 
#include "newhash.h"
#include "extfunc.h" 
#include "extvab.h" 

static int Ncounter;
static int allGaps;

// for multi threads
static int scafBufSize = 100;
static STACK **ctgStackBuffer;
static int scafCounter;
static int scafInBuf;

static void convertIndex()
{
	int *length_array = (int *)ckalloc((num_ctg+1)*sizeof(int));
	unsigned int i;
	for(i=1;i<=num_ctg;i++)
		length_array[i] = 0;

	for(i=1;i<=num_ctg;i++){
		if(index_array[i]>0)
			length_array[index_array[i]] = i;
	}
	for(i=1;i<=num_ctg;i++)
		index_array[i] = length_array[i];  //contig i with new index: index_array[i]
	free((void *)length_array);

}

static void reverseStack(STACK *dStack,STACK *sStack)
{
	CTGinSCAF *actg,*ctgPt;
	emptyStack(dStack);

	while((actg=(CTGinSCAF *)stackPop(sStack))!=NULL){
		ctgPt = (CTGinSCAF *)stackPush(dStack);
		ctgPt->ctgID = actg->ctgID;
		ctgPt->start = actg->start;
		ctgPt->end = actg->end;
	}
	stackBackup(dStack);
}

static void initStackBuf(STACK **ctgStackBuffer,int scafBufSize)
{
	int i;
	for(i=0;i<scafBufSize;i++)
		ctgStackBuffer[i] = (STACK *)createStack(100,sizeof(CTGinSCAF));

}
static void freeStackBuf(STACK **ctgStackBuffer,int scafBufSize)
{
	int i;
	for(i=0;i<scafBufSize;i++)
		freeStack(ctgStackBuffer[i]);
}

static void mapCtg2Scaf(int scafInBuf)
{
	int i,scafID;
	CTGinSCAF *actg;
	STACK *ctgsStack;
	unsigned int ctg,bal_ctg;
	
	for(i=0;i<scafInBuf;i++){
		scafID = scafCounter+i+1;
		ctgsStack = ctgStackBuffer[i];
		while((actg=stackPop(ctgsStack))!=NULL){
			ctg = actg->ctgID;
			bal_ctg = getTwinCtg(ctg);
			
			if(contig_array[ctg].from_vt!=0){
				contig_array[ctg].multi = 1;
				contig_array[bal_ctg].multi = 1;
				continue;
			}

			contig_array[ctg].from_vt = scafID;
			contig_array[ctg].to_vt = actg->start;
			contig_array[ctg].flag = 0;   //ctg and scaf on the same strand
			contig_array[bal_ctg].from_vt = scafID;
			contig_array[bal_ctg].to_vt = actg->start;
			contig_array[bal_ctg].flag = 1;
		}
	}

}

static void locateContigOnscaff(char *graphfile)
{

	FILE *fp;
	char line[1024];
	CTGinSCAF *actg;
	STACK *ctgStack,*aStack;
	int index=0,counter,overallLen;
	int starter,prev_start,gapN,scafLen = 0;
	unsigned int ctg,prev_ctg=0;

	for(ctg=1;ctg<=num_ctg;ctg++){
		contig_array[ctg].from_vt = 0;
		contig_array[ctg].multi = 0;
	}

	ctgStack = (STACK *)createStack(1000,sizeof(CTGinSCAF));

	sprintf(line, "%s.scaf_gap", graphfile);
	fp = ckopen(line, "r");

	ctgStackBuffer = (STACK **)ckalloc(scafBufSize*sizeof(STACK *));
	initStackBuf(ctgStackBuffer,scafBufSize);


	Ncounter = scafCounter = scafInBuf = allGaps = 0;
	while(fgets(line,sizeof(line),fp)!=NULL){
		if(line[0]=='>'){
			if(index){
				aStack = ctgStackBuffer[scafInBuf++];
				reverseStack(aStack,ctgStack);
				if(scafInBuf==scafBufSize){
					mapCtg2Scaf(scafInBuf);
					scafCounter += scafInBuf;
					scafInBuf = 0;
				}
				//if(index%1000==0)
					//printf("Processed %d scaffolds\n",index);
			}
			//read next scaff	
			scafLen = prev_ctg = 0;
			emptyStack(ctgStack);
			sscanf(line+9,"%d %d %d",&index,&counter,&overallLen);
			fprintf(stderr,">%d\n",index);
			continue;
		}
		if(line[0]=='G'){  // gap appears
			continue;
		}
		if(line[0]>='0'&&line[0]<='9'){   // a contig line
			sscanf(line,"%d %d",&ctg,&starter);
			actg = (CTGinSCAF *)stackPush(ctgStack);	
			actg->ctgID = ctg;
			if(!prev_ctg){
				actg->start = scafLen;
				actg->end = actg->start + overlaplen + contig_array[ctg].length - 1;
			}else{
				gapN = starter - prev_start-(int)contig_array[prev_ctg].length;
				gapN = gapN < 1 ? 1:gapN;
				actg->start = scafLen + gapN;
				actg->end = actg->start + contig_array[ctg].length - 1;
			}
			fprintf(stderr,"%d\t%d\n",actg->start,actg->end);
			scafLen = actg->end+1;
			prev_ctg = ctg;
			prev_start = starter;
		}
	}
	if(index){
		aStack = ctgStackBuffer[scafInBuf++];
		reverseStack(aStack,ctgStack);
		mapCtg2Scaf(scafInBuf);
	}
	gapN = 0;
	for(ctg=1;ctg<=num_ctg;ctg++){
		if(contig_array[ctg].from_vt==0||contig_array[ctg].multi==1)
			continue;
		gapN++;
	}
	//printf("\nDone with %d scaffolds, %d contigs in Scaffolld\n",index,gapN);
	fclose(fp);
	freeStack(ctgStack);
	freeStackBuf(ctgStackBuffer,scafBufSize);
	free((void*)ctgStackBuffer);
}

static boolean contigElligible(unsigned int contigno)
{
	unsigned int ctg = index_array[contigno];
	if(contig_array[ctg].from_vt==0||contig_array[ctg].multi==1)
		return 0;	
	else
		return 1;	

}
static void output1read(FILE *fo,long long readno,unsigned int contigno,int pos)
{
	
	unsigned int ctg = index_array[contigno];
	int posOnScaf;
	char orien;
	pos = pos < 0 ? 0:pos;
	if(contig_array[ctg].flag==0){
		posOnScaf = contig_array[ctg].to_vt + pos - overlaplen;
		orien = '+';
	}else{
		posOnScaf = contig_array[ctg].to_vt + contig_array[ctg].length - pos;
		orien = '-';
	}
	/*
	if(readno==676)
		printf("Read %lld in region from %d, extend %d, pos %d, orien %c\n",
			readno,contig_array[ctg].to_vt,contig_array[ctg].length,posOnScaf,orien);
	*/
	fprintf(fo,"%lld\t%d\t%d\t%c\n",readno,contig_array[ctg].from_vt,posOnScaf,orien);
}

void locateReadOnScaf(char *graphfile)
{
	char name[1024],line[1024];
	FILE *fp,*fo;
	long long readno,counter=0,pre_readno=0;
	unsigned int contigno,pre_contigno = 0;
	int pre_pos = 0, pos;

	locateContigOnscaff(graphfile);

	sprintf(name,"%s.readOnContig",graphfile);
	fp = ckopen(name,"r");
	sprintf(name,"%s.readOnScaf",graphfile);
	fo = ckopen(name,"w");

	if(!orig2new){
		convertIndex();
		orig2new = 1;
	}
	fgets(line,1024,fp);
	while(fgets(line,1024,fp)!=NULL){
		sscanf(line,"%lld %d %d",&readno,&contigno,&pos);
		if((readno%2==0)&&(pre_readno==readno-1) // they are a pair of reads
			&&contigElligible(pre_contigno)&&contigElligible(contigno)){
			output1read(fo,pre_readno,pre_contigno,pre_pos);
			output1read(fo,readno,contigno,pos);
			counter++;
		}
		pre_readno = readno;
		pre_contigno = contigno;
		pre_pos = pos;
	}
	printf("%lld pairs on contig\n",counter);
	fclose(fp);
	fclose(fo);
}
